% PSO Optimization Script
clear;
clc;

% 示例数据集
data = {
    56769, '2015-04-17', '17:17:30', '6号线', '金桥路', '17:24:37', '6号线', '博兴路', 3.0;
    56769, '2015-04-17', '19:59:57', '6号线', '博兴路', '20:04:07', '6号线', '金桥路', 3.0;
    %... （其他数据行）
    427031, '2015-04-17', '10:51:19', '6号线', '巨峰路', '11:03:06', '6号线', '金桥路', 2.7
};

% 将时间转换为秒数
for i = 1:size(data, 1)
    data{i, 3} = time_to_seconds(data{i, 3});
    data{i, 6} = time_to_seconds(data{i, 6});
end

% 转换数据为矩阵
data = cell2mat(data(:, [1, 3, 6, 9]));
travel_time = data(:, 3) - data(:, 2);  % 计算出行时间

% PSO算法参数
num_particles = 50;
max_iter = 80;
dim = 3;
lb = -10 * ones(1, dim);
ub = 10 * ones(1, dim);
w_max = 0.9;
w_min = 0.4;
c1 = 2;
c2 = 2;

% 初始化粒子位置和速度
pos = lb + (ub - lb) .* rand(num_particles, dim);
vel = zeros(num_particles, dim);

% 初始化个体最优位置和全局最优位置
pbest = pos;
pbest_val = zeros(num_particles, 1);
params = struct('alpha', 1, 'beta', 1, 'gamma', 1, 'C_b', 1, 'H_min', 18000, 'H_max', 79200, 'C_k', 50); % 示例参数
for i = 1:num_particles
    [pbest_val(i), ~, ~, ~] = objective_function(pbest(i, :), params, travel_time);
end
[~, gbest_idx] = min(pbest_val);
gbest = pbest(gbest_idx, :);
gbest_val = pbest_val(gbest_idx);

% 存储每次迭代的最优值
best_values = zeros(max_iter, 1);
best_Cx = zeros(max_iter, 1);
best_Cy = zeros(max_iter, 1);
best_Cv = zeros(max_iter, 1);

% PSO主循环
for iter = 1:max_iter
    w = w_max - (w_max - w_min) * (iter / max_iter);
    for i = 1:num_particles
        vel(i, :) = w * vel(i, :) + c1 * rand() * (pbest(i, :) - pos(i, :)) + c2 * rand() * (gbest - pos(i, :));
        pos(i, :) = pos(i, :) + vel(i, :);
        pos(i, :) = max(min(pos(i, :), ub), lb);
        [current_val, Cx, Cy, Cv] = objective_function(pos(i, :), params, travel_time);
        if current_val < pbest_val(i)
            pbest(i, :) = pos(i, :);
           pbest_val(i) = current_val;
        end
        if current_val < gbest_val
            gbest = pos(i, :);
            gbest_val = current_val;
            best_Cx(iter) = Cx;
            best_Cy(iter) = Cy;
            best_Cv(iter) = Cv;
        end
    end
    best_values(iter) = gbest_val;
    fprintf('Iter %d: Best Value = %.4f, Cx = %.4f, Cy = %.4f, Cv = %.4f\n', iter, gbest_val, best_Cx(iter), best_Cy(iter), best_Cv(iter));
end

% 结果展示和图形化
disp('Optimization completed.');
disp(['Best Value: ', num2str(gbest_val)]);
disp(['Best Position: ', num2str(gbest)]);
figure;
subplot(4, 1, 1);
plot(1:max_iter, best_values, '-o');
xlabel('Iteration');
ylabel('Best Value');
title('PSO Optimization Progress');
subplot(4, 1, 2);
plot(1:max_iter, best_Cx, '-o');
xlabel('Iteration');
ylabel('C_x');
title('C_x Progress');
subplot(4, 1, 3);
plot(1:max_iter, best_Cy, '-o');
xlabel('Iteration');
ylabel('C_y');
title('C_y Progress');
subplot(4, 1, 4);
plot(1:max_iter, best_Cv, '-o');
xlabel('Iteration');
ylabel('C_v');
title('C_v Progress');
grid on;

% 时间转换函数
function seconds = time_to_seconds(time_str)
    t = datetime(time_str, 'InputFormat', 'HH:mm:ss');
    seconds = hour(t) * 3600 + minute(t) * 60 + second(t);
end

% 目标函数定义
function [Z, Cx, Cy, Cv] = objective_function(X, params, travel_time)
    alpha = params.alpha;
    beta = params.beta;
    gamma = params.gamma;
    C_b = params.C_b;
    H_min = params.H_min;
    H_max = params.H_max;
    C_k = params.C_k;
    
    % 从X向量中提取具体变量
    A_nk = X(1:end-2);
    T_Ai_k = X(end-1);
    T_Li_k = X(end);
    
    % 计算各部分成本
    C_d = abs(T_Ai_k - T_Li_k);
    Cx = C_b + C_d;
    T_P_wait = abs(T_Ai_k - T_Li_k);
    Cy = C_b + T_P_wait;
    Cv = T_P_wait;
    
    % 添加约束条件
    penalty = 0;
    
    % 约束1: 列车计划约束
    if T_Ai_k < H_min || T_Li_k > H_max
        penalty = penalty + 1e6;
    end
    
    % 约束2: 列车容量限制约束
    if sum(A_nk) > C_k
        penalty = penalty + 1e6;
    end
    
    % 约束3: 乘客服务约束
    if sum(A_nk) ~= 1
        penalty = penalty + 1e6;
    end
    
    % 目标函数
    Z = alpha * sum(Cx .* travel_time) + beta * sum(Cy .* travel_time) + gamma * sum(Cv .* travel_time) + penalty;
end
